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Abstract 

Phase oscillators are a common starting point for the reduced description of many single neuron models that 
exhibit a strongly attracting limit cycle. The framework for analysing such models in response to weak 
perturbations is now particularly well advanced, and has allowed for the development of a theory of weakly 
connected neural networks. However, the strong-attraction assumption may well not be the natural one for 
many neural oscillator models. For example, the popular conductance based Morrls-Lecar model is known to 
respond to periodic pulsatile stimulation in a chaotic fashion that cannot be adequately described with a phase 
reduction. In this paper, we generalise the phase description that allows one to track the evolution of distance 
from the cycle as well as phase on cycle. We use a classical technique from the theory of ordinary differential 
equations that makes use of a moving coordinate system to analyse periodic orbits. The subsequent 
phase-amplitude description is shown to be very well suited to understanding the response of the oscillator to 
external stimuli (which are not necessarily weak). We consider a number of examples of neural oscillator models, 
ranging from planar through to high dimensional models, to illustrate the effectiveness of this approach in 
providing an improvement over the standard phase-reduction technique. As an explicit application of this 
phase-amplitude framework, we consider in some detail the response of a generic planar model where the 
strong-attraction assumption does not hold, and examine the response of the system to periodic pulsatile 
forcing. In addition, we explore how the presence of dynamical shear can lead to a chaotic response. 



1 



Keywords 

phase-amplitude, oscillator, chaos, non-weak coupling 
1 Introduction 

One only has to look at the plethora of papers and books on the topic of phase oscillators in mathematical 
neuroscience to see the enormous impact that this tool from dynamical systems theory has had on the way 
we think about describing neurons and neural networks. Much of this work has its roots in the theory of 
ordinary differential equations (ODEs) and has been promoted for many years in the work of Winfree [l], 
Guckenheimcr [2], Holmes Kopell Ermcntrout 5| and Izhikevich |6 to name but a few. For a recent 
survey we refer the reader to the book by Ermcntrout and Terman |7 . At heart, the classic phase reduction 
approach recognises that if a high dimensional nonlinear dynamical system (as a model of a neuron) 
exhibits a stable limit cycle attractor then trajectories near the cycle can be projected onto the cycle. 
A natural phase variable is simply the time along the cycle (from some arbitrary origin) relative to the 
period of oscillation. The notion of phase can even be extended off the cycle using the concept of 
isochrons |1 . They provide global information about the 'latent phase', namely the phase that will be 
asymptotically returned to for a trajectory with initial data within the basin of attraction of an 
exponentially stable periodic orbit. More technically, isochrons can be viewed as the leaves of the invariant 
foliation of the stable manifold of a periodic orbit |8 . In rotating frame coordinates given by phase and the 
leaf of the isochron foliation the system has a skew-product structure, i.e. the equation of the phase 
decouples. However, it is a major challenge to find the isochron foliation, and since it relies on the 
knowledge of the limit cycle it can only be found in special cases or numerically. There are now a number 
of complementary techniques that tackle this latter challenge, and in particular we refer the reader to work 
of Guillamon and Huguet [oj (using Lie symmetries) and Osinga and Moehlis [lO] (exploiting numerical 
continuation). More recent work by Mauroy and Mezic [ll] is especially appealing as it uses a simple 
forward integration algorithm, as illustrated in Fig. [l] for a Stuart-Landau oscillator. However, it is more 
common to side-step the need for constructing global isochrons by restricting attention to a small 
neighbourhood of the limit cycle, where dynamics can simply be recast in the reduced form = 1, where 9 
is the phase around a cycle. This reduction to a phase description gives a nice simple dynamical system, 
albeit one that cannot describe evolution of trajectories in phase-space that are far away from the limit 
cycle. However, the phase reduction formalism is useful in quantifying how a system (on or close to a 
cycle) responds to weak forcing, via the construction of the infinitesimal phase response curve (iPRC). For 
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a given high dimensional conductance based model this can be solved for numerically, though for some 
normal form descriptions closed form solutions are also known 12 . The iPRC at a point on cycle is equal 
to the gradient of the (isochronal) phase at that point. This approach forms the basis for constructing 
models of weakly interacting oscillators, where the external forcing is pictured as a function of the phase of 
a firing neuron. This has led to a great deal of work on phase-locking and central pattern generation in 



neural circuitry (see for example 13 ). Note that the work in |9j goes beyond the notion of iPRC and 
introduces infinitesimal phase response surfaces (allowing evaluation of phase advancement even when the 
stimulus is off cycle), and see also the work in (14| on non- infinitesimal PRCs. 




Figure 1: Isochrons of a Stuart-Landau oscillator model: x = Xx/2 — (Ac/2 -|- uj)y — X{x^ + y'^){x — cy)/2, 
y = (Ac/2 -I- uj)x + Ay/2 — X{x^ + y'^){cx + y)/2. The black curve represents the periodic orbit of the 
system, which is simply the unit circle for this model. The blue curves are the isochrons obtained using 
the (forward) approach of Mauroy and Mezic [iT]. The green dots are analytically obtained isochronal 
points [15]. Parameter values are A = 2.0, c = 1.0 and w = 1.0. 



The assumption that phase alone is enough to capture the essentials of neural response is one made more 
for mathematical convenience than being physiologically motivated. Indeed for the popular type I 
Morris-Lecar (ML) firing model with standard parameters, direct numerical simulations with pulsatile 



forcing show responses that cannot be explained solely with a phase model 16 . The failure of a phase 



description is in itself no surprise and underlies why the community emphasises the use of the word weakly 
in the phrase "weakly connected neural networks" . Indeed, there are a number of potential pitfalls when 
applying phase reduction techniques to a system that is not in a weakly forced regime. The typical 
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construction of the phase response curve uses only hnear information about the isochrons and nonhnear 
effects win come into play the further we move away from the limit cycle. This problem can be diminished 
by taking higher order approximations to the isochrons and using this information in the construction of a 
higher order PRC 17]. Even using perfect information about isochrons, the phase reduction still assumes 
persistence of the limit-cycle and instantaneous relaxation back to cycle. However, the presence of nearby 
invariant phase-space structures such as (unstable) fixed point and invariant manifolds may result in 
trajectories spending long periods of time away from the limit cycle. Moreover, strong forcing will 
necessarily take one away from the neighbourhood of a cycle where a phase description is expected to hold. 
Thus developing a reduced description which captures some notion of distance from cycle is a key 
component of any theory of forced limit cycle oscillators. The development of phase-amplitude models that 
better characterise the response of popular high dimensional single neuron models is precisely the topic of 
this paper. Given that it is a major challenge to construct an isochronal foliation we use non-isochronal 
phase-amplitude coordinates as a practical method for obtaining a more accurate description of neural 



systems. Recently Medvedev 18 has used this approach to understand in more detail the synchronisation 
of linearly coupled stochastic limit cycle oscillators. 

In section [2| we consider a general coordinate transformation which recasts the dynamics of a system in 
terms of phase-amplitude coordinates. This approach is directly taken from the classical theory for 



analysing periodic orbits of ODEs, originally considered for planar systems in 19 , and for general systems 



20 . We advocate it here as one way to move beyond a purely phase-centric perspective. We illustrate 



the transformation by applying it to a range of popular neuron models. In section [3j we consider how 
inputs to the neuron are transformed under these coordinate transformations and derive the evolution 
equations for the forced phase-amplitude system. This reduces to the standard phase description in the 
appropriate limit. Importantly, we show that the behaviour of the phase- amplitude system is much more 
able to capture that of the original single neuron model from which it is derived. Focusing on pulsatile 



forcing we explore the conditions for neural oscillator models to exhibit shear induced chaos 16 . Finally in 
section |4j we discuss the relevance of this work to developing a theory of network dynamics that can 
improve upon the standard weak coupling approach. 



2 Phase-amplitude coordinates 

Throughout this paper, we study the dynamics prescribed by the system x — f{x), x € M", with solutions 
X = x{t) that satisfy a;(0) — xq ^ M". We will assume that the system admits an attracting hyperbolic 
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periodic orbit (namely one zero Floquet exponent and the others having negative real part), with period A, 
such that u{t) = u{t + A). A phase 6 e [0, A) is naturally defined from 9{u{t)) = t mod A. It has long 
been known in the dynamical systems community how to construct a coordinate system based on this 



notion of phase as well as a distance from cycle, see 20 for a discussion. In fact, Ermentrout and 



Kopell 21 made good use of this approach to derive the phase-interaction function for networks of weakly 
connected limit-cycle oscillators in the limit of infinitely fast attraction to cycle. However, this assumption 
is particularly extreme and unlikely to hold for a broad class of single neuron models. Thus, it is 
interesting to return to the full phase-amplitude description. In essence, the transformation to these 
coordinates involves setting up a moving orthonormal system around the limit cycle. One axis of this 
system is chosen to be in the direction of the tangent vector along the orbit, and the remaining are chosen 
to be orthogonal. We introduce the normalised tangent vector ^ as 



m = 



Au 



Au 
A0 



(1) 



The remaining coordinate axes are conveniently grouped together as the columns of an n x (n — 1) matrix 
Q. In this case we can write an arbitrary point x as 



(2) 



Here, \p\ represents the Euclidean distance from the limit cycle. A caricature, in , of the coordinate 
system along an orbit segment is shown in Fig. [2] Through the use of the variable p, we can consider 
points away from the periodic orbit. Rather than being isochronal, lines of constant 9 are simply straight 
lines that emanating from a point on the orbit in the direction of the normal. The technical details of 
specifying the orthonormal coordinates forming C, are discussed in Appendix [A] 
Upon projecting the dynamics onto the moving orthonormal system, we obtain the dynamics of the 
transformed system: 



= l + p^A{9)p + h{9,p). 



where 



hi 

/2( 



, p) = ^h^i9, p)^p + h^{9, p) [J{u + Cp) - f{u)] 



h{9,p) = 



Au 
A9 



t Ai9)=e 



A9 



D/C 



(3) 

(4) 
(5) 

(6) 



5 




Figure 2: Demonstration of the moving orthonormal coordinate system along an orbit segment. As t evolves 
from to to ti, the coordinates vary smoothly. In this planar example, C always points to the outside of the 
orbit. 



and Vf is the Jacobian of the vector field /, evaluated along the periodic orbit u. The derivation of this 
system may be found in appendix [A| It is straight-forward to show that fi{0, p) — > as \p\ — > 0, f2iS, 0) = 
and that 8/2(0, 0)/dp — 0. In the above, /i captures the shear present in the system, that is, whether the 
speed of 9 increases or decreases dependent on the distance from cycle. A precise definition for shear is 



given in 22 . Additionally, A{9) describes the ^-dependent rate of attraction or repulsion from cycle. 
It is pertinent to consider where this coordinate transformation breaks down, that is, where the 
determinant of the Jacobian of the transformation K ~ det[dx/d6 dx/dp] vanishes. This never vanishes 
on-cycle (where p = 0), but may do so for some \p\ = k > 0. This sets an upper bound on how far away 
from the limit cycle we can describe the system using these phase-amplitude coordinates. In Fig. [3] we plot 
the curve along which the transformation breaks down for the ML model. We observe that, for some values 
of 9, k is relatively smaller. The breakdown occurs where lines of constant 9 cross and thus the 
transformation ceases to be invertible, and these values of correspond to points along which the orbit has 
high curvature. We note that this issue is less problematical in higher dimensional models. 
If we now consider the driven system 

x^ f{x)+eg{x,t), xeM", (7) 

where e is not necessarily small, we may apply the same transformation as above to obtain the dynamics in 
{9, p) coordinates, where e [0, A) and p £ M"^^, as 

9 = 1 + M9,p)+ eh^{9, p)g{u{9) + C(0)p, i), (8) 
p = A{9)p + f2{9, p) + e(^B{9, p)g{u{9) + CWp, t), (9) 
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Figure 3: This figure shows the determinant of the phase-amphtude transformation for the ML modeL 
Colours indicate the value of K. The red contour indicates where K = 0, and thus where the coordinate 
transformation breaks down. Note how all of the values for which this occurs have p < 0. Parameter values 



as in Appendix B.l 



with 

BiO,p)^ln~'^ph^{0.p), (10) 

and I„ is the n x n identity matrix. Here, h and B describe the effect in terms of 6 and p that the 
perturbations have. Details of the derivation are given in Appendix [X] For planar models, S = I2. To 
demonstrate the application of the above coordinate transformation, we now consider some popular single 
neuron models. 

2.1 A 2D conductance based model 

The ML model was originally developed to describe the voltage dynamics of barnacle giant muscle 

fiber j23|, and is now a popular modelling choice in computational neuroscience |7l. It is written as a pair 

of coupled nonlinear ODEs of the form 

Cv = I{t) ~ gi{v ~ Vl) - 5kw(w - Vk) - 5caTOoo(w)(w - Wca), W = (f>{Wao{v) - w)/t.u,{v). (11) 

Here, v is the membrane voltage, whilst w is a gating variable^ describing the fraction of membrane ion 
channels that are open at any time. The first equation expresses Kirchoff 's current law across the cell 
membrane, with I(t) representing a stimulus in the form of an injected current. The detailed form of the 
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Figure 4: Typical trajectory of the ML model of the transformed system. Left: Time evolution of 9 and p. 
Right: Trajectory plotted in the (v, w) phase plane. We see that when p has a local maximum, the evolution 
of 6 slows down, corresponding to where trajectories pass near to the saddle node. Parameter values as in 
Appendix |B.1[ 



model is completed in Appendix |B.1[ The ML model has a very rich bifurcation structure. Roughly 
speaking, by varying a constant current I{t) = Iq , one observes, in different parameter regions, dynamical 
regimes corresponding to sinks, limit cycles, and Hopf, saddle-node and homoclinic bifurcations, as well as 
combinations of the above. These scenarios are discussed in detail in |7] and 24 . 



As the ML model is planar, p is a scalar, as are the functions A and /i.2. This allows us to use the moving 
coordinate system to clearly visualise parts of phase space where trajectories are attracted towards the 
limit cycle, and parts in which they move away from it, as illustrated in Fig. ^ The functions /1.2 and A, 
evaluated at p — —0.1 are shown in Fig. [5] The evolution of 6 is mostly constant, however we clearly 
observe portions of the trajectories where this is slowed, along which, p ~ 0. In fact, this corresponds to 
where trajectories pass near to the saddle node, and the dynamics stall. This occurs around 6 = 17, and in 
Fig. [5] we see that both A{9) and fi{0, p) are indeed close to 0. The reduced velocities of trajectories here 
highlights the importance of considering other phase space structures in forced systems, the details of 
which arc missed in standard phase only models. Forcing in the presence of such structures may give rise 
to complex and even chaotic behaviours, as we shall see in section [3j 

In the next example, we show how the same ideas go across to higher dimensional models. 
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Figure 5: /i,/2 and A for the ML model, evaluated at p = —0.1. We clearly see the difference in the 
order of magnitude between /i and /2 for small p. Note that, although the average of A over one period 
is negative, it is positive for a nontrivial interval of 9. This corresponds to movement close to the stable 
manifold of the saddle node. Parameter values as in Appendix |B.1| 



2.2 A 4D conductance based model 



The Connor-Stevens (CS) model 25 is built upon the Hodgkin-Huxley formalism and comprises a fast 
Na^ current, a delayed current, a leak current and a transient current, termed the A-current. The 
full CS model consists of 6 equations: the membrane potential, the original Hodgkin-Huxley gating 
variables and an activating and inactivating gating variable for the A-current. Using the method of 



equivalent potentials 26 , we may reduce the dimensionality of the system, to include only 4 variables. The 
reduced system is: 

Xoo{v)-X 



Cv = —F{v,u,a,b) + I, u — G{v,u), X = 



Tx{v) 



Xe{a,b}, 



(12) 
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where 

F{v, u, a, b) = gKnf^{u){v - vk) + 5Na^oo('«)?nL(")(" - "^^Na) + ghKv - wi) + ge,a^b{v - v^). (13) 

The details of the reduced CS model are completed in Appendix |B.2| The solutions to the reduced CS 
model under the coordinate transformation may be seen in Fig.|6l whilst, in Fig. [71 we show how this 
solution looks in the original coordinates. As for the ML model, 9 evolves approximately constantly 
throughout, though this evolution is sped up close to 6* = A. The trajectories of the vector p arc more 
complicated, but note that there is regularity in the pattern exhibited, and that this occurs with 
approximately the same period as the underlying limit cycle. The damping of the amplitude of oscillations 
in p over successive periods represents the overall attraction to the limit cycle, whilst the regular behaviour 
of p represents the specific relaxation to cycle as shown in Fig. [Tj Additional file 1 shows a movie of the 
trajectory in (w, u, h) coordinates with the moving orthonormal system superimposed, as well as the 
solution in p for comparison. 



3 Pulsatile forcing of phase-amplitude oscillators 

We now consider a system with time-dependent forcing, given by ([t]) with 

5(x,0 = ^(^(i-nr),0,...,0)^, 



(14) 



where 5 is the Dirac (5-function. This describes T-periodic kicks to the voltage variable. Even such a simple 



forcing paradigm can give rise to rich dynamics 16 . For the periodically kicked ML model, shear forces 



can lead to chaotic dynamics as folds and horseshoes accumulate under the forcing. This means that the 
response of the neuron is extremely sensitive to the initial phase when the kicks occur. In terms of neural 



response, this means that the neuron is unreliable 27 



The behaviour of oscillators under such periodic pulsatile forcing is the subject of a number of studies, see, 
e-g-, 



27 -30 . Of particular relevance here is 27 , in which a qualitative reasoning of the mechanisms that 



bring about shear in such models is supplemented by direct numerical simulations to detect the presence of 
chaotic solutions. For the ML model in a parameter region close to the homoclinic regime, kicks can cause 



trajectories to pass near the saddle-node, and folds may occur as a result 16 



We here would like to compare full planar neural models to the simple model, studied in 27 



= l + cr/9, p = ~\p + eP{e)^5{t-nT). 



(15) 
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Figure 6: Solution of the transformed CS system. Top: Time evolution of 6. Bottom: Time evolution of p 
coordinates. Upon transforming these coordinates back to the original ones, we arrive at Fig. [7j Parameter 
values given in Appendix |B.2| In this parameter regime the model exhibits type I firing dynamics. 



This system exhibits dynamical shear which, under certain conditions, can lead to chaotic dynamics. The 
shear parameter cr dictates how much trajectories are 'sped up' or 'slowed down' dependent on their 
distance from the limit cycle, whilst A is the rate of attraction back to the limit cycle, which is independent 
of 9. Supposing that the function P is smooth but non-constant, trajectories will be taken a variable 
distance from the cycle upon the application of the kick. When kicks are repeated, this geometric 
mechanism can lead to repeated stretching and folding of phase space. It is clear that the larger a is in 



(15), the more shear is present, and the more likely we are to observe the folding effect. In a similar way. 



smaller values of A mean that the shear has longer to act upon trajectories and again result in a greater 
likelihood of chaos. Finally, to observe chaotic response, we must ensure that the shear forces have 
sufficient time to act, meaning that T, the time between kicks must not be too small. 
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Figure 7: Transformed trajectory in {v,u,b) space of the phase-amplitude description of the reduced CS 
model. The dotted black line is the phase amplitude solution transformed in the original coordinates, whilst 
the coloured orbit is the underlying periodic orbit, where the colour corresponds to the phase along the orbit. 
The solution of the phase-amplitude description of the model in {6, p) coordinates is shown in Fig. p] 



This stretching and folding action can clearly lead to the formation of Smale horseshoes, which are well 
known to lead to a type of chaotic behaviour. However, horseshoes may coexist with sinks, meaning the 
resulting chaotic dynamics would be transient. Wang and Young proved that under appropriate conditions, 
there is a set of T of positive Lebesgue measure for which the system experiences a stronger form of 
sustained, chaotic behaviour, characterised by the existence of a positive Lyapunov exponent for almost all 



initial conditions and the existence of a "strange attractor"; see, e.g., 28 -30 



By comparing with the phase-amplitude dynamics described by equations (|8])-(|9|, we see that the model of 
shear considered in ( |15[ ) is a proxy for a more general system, with fi{9,p) — > ap, A{6) — > —A and 
h{e,p) -> 0, and C{d) P{d)- 

To gain a deeper insight into the phenomenon of shear induced chaos, it is pertinent to study the isochrons 
of the limit cycle for the linear model (15), where the isochrons are simply lines with slope —\/a. In Fig. [sj 
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Figure 8: Stretch-and-fold action of a kick followed by relaxation in the presence of shear. The thin black 
lines are the isochrons of the system, which in the case of the linear model ( 15 ), are simply straight lines with 
slope The thin gray line at p 

with strength e 



represents the limit cycle, which is kicked, aXt — by P{9) = sin(6') 
1 to the solid curve. After this, the orbits are allowed to evolve under the flow generated 
by the continuous part of the system. The dashed and dotted curves represent the image of the kicked solid 
curve under this flow, at times ti and ^2 respectively. The green marker shows how one point, a;(to) evolves 
under the flow, first to x{ti) and then to x{t2), following the isochron as it relaxes back to the limit cycle. 
The effect of the shear forces and the subsequent folding, caricatured by the blue arrows can clearly be seen. 



we depict the isochrons and stretch and fold action of shear. The bold thin gray line at p = is kicked, at 
t ^ to, to the bold solid curve, where P{9) ~ sin(6'), as studied in [l6] and this curve is allowed to evolve 
under the dynamics with no further kicks through the dashed curve at < = ii and ultimately to the dotted 
curve at i = ^2 , which may be considered as evolutions of the solid curve by integer multiples of the natural 
period of the oscillator. Every point of the dotted curve traverses the isochron it is on at until <2- The 
green marker shows an example of one such point evolving along its associated isochron. The folding effect 
of this is clear in the figure, and further illustrated in the video in Additional file 2. 
This simple model with a harmonic form for P(&) provides insight into how strange attractors can be 
formed. Kicks along the isochrons or ones that map isochrons to one another will not produce strange 
attractors, but merely phase-shifts. What causes the stretching and folding is the variation in how far 



points are moved as measured in the direction transverse to the isochrons. For the linear system (15 1 
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variation in this sense is generated by any non-constant P{0); the larger the ratio ere/ A, the larger the 
variation (see [16] for a recent discussion). 

The formation of chaos in the ML model is shown in Fig. [9| Here, we plot the response to periodic 



pulsatile forcing, given by (14), in the {v,w) coordinate system. This clearly illustrates a folding of phase 



space around the limit-cycle, and is further portrayed in the video in Additional file 3. We now show how 
this can be understood using phase-amplitude coordinates. 
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Figure 9: Shear induced folding in the ML model, with parameters as in Fig. [5| The red curve in all panels 
represents the limit cycle of the unperturbed system, whilst the green dotted line represents the stable 
manifold of the saddle node, indicated by the orange marker. We begin distributing points along the limit 
cycle and then apply an instantaneous kick taking v v + e, where e = —2.0, leaving w unchanged. This 
essentially moves all phase points to the left, to the blue curve. The successive panels show the image of this 
set of points after letting points evolve freely under the system defined by the ML equations, and then apply 
the kick again. The curves shown are the images of the initial phase points just after each kick, as indicated 
in the figure. We can clearly observe the shear induced folding. Parameter values as in Appendix |B.1| 



Compared to the phenomenological system (15), models written in phase-amplitude coordinates as (^8|9 



14 



have two main differences. The intrinsic dynamics (without kicks) are nonhnear, and the kick terms appear 
in both equations for 9 and p (not just p). Simulations of ([8])-([9| for both the FHN and ML models, with 
e — 0.1, show that the replacement of fi{9,p) by ap, dropping f2{0,p) (which is quadratic in p), and 
setting A{9) ~ —A does not lead to any significant qualitative change in behaviour (for a wide range of 
cr, A > 0). We therefore conclude that, at least when the kick amplitude e is not too large, it is more 
important to focus on the form of the forcing in the phase-amplitude coordinates. In what follows, we are 
interested in discovering the effects of different functional forms of the forcing term multiplying the 
(5-function, keeping other factors fixed. As examples, we choose those forcing terms given by transforming 
the FHN and the ML models into phase-amplitude coordinates. To find these functions, we first find the 



attracting limit cycle solution in the ML model (111 and FHN model (52 1 using a periodic boundary value 
problem solver and set up the orthonormal coordinate system around this limit cycle. Once the coordinate 
system is established, we evaluate the functions h{9,p) and B{9,p) (that appear in equations (Is]) and (|9|). 



For planar systems we have simply that B{9,p) = I2. Using the forcing term (14), we are only considering 
perturbations to the voltage component of our system and thus only the first component of h and B will 
make a nontrivial contribution to the dynamics. We define Pi as the first component of h and P2 as the first 
component of We wish to force each system at the same ratio of the natural frequency of the underlying 
periodic orbit. To ease comparison between the system with the ML forcing terms and the FHN forcing 
terms, we rescale 9 1-^ 9/ A so that 9 E [0, 1) in what follows. Implementing the above choices leads to 

9 ^l + ap + sPi{9,p)J2^it-nT), p^-\p + eP2{9)^5{t-nT). (16) 
It is important to emphasise that Pi ^2 are determined by the underlying single neuron model (unlike in the 



toy model (15l). As emphasised in 31 , one must take care in the treatment of the state dependent 'jumps' 



caused by the (5-functions in ( 16 ) to accommodate the discontinuous nature of 9 and p at the time of the 



kick. To solve (16), we approximate 5{t) with a normalised square pulse (^^-(t) of the form 



[0 i<0, 

5r{t)=\l/T 0<t<T, (17) 
[0 t>T, 

where t ^ 1. This means that for (n — 1)T + t < t < nT, n E Z, the dynamics are governed by the linear 
system {9,p) = (1 + <jp, — Ap). This can be integrated to obtain the state of the system just before the 



15 



arrival of the next kick, 



3'n,Pri) = {S{nT),p{nT)), in the form 

= \e{t) + nT~t + jp{t) (i - e-^("^-*)) 



mod 1, 



(18) 
(19) 



In the interval nT < t < nT + r and using (17 1 we now need to solve the system of nonlinear ODEs 



(20) 



Rescaling time as t = nT + rs, with < s < 1, and writing the solution (9, p) as a regular perturbation 
expansion in powers of r as {9{s), p{s)) — {do{s) + t9i(s), po{s) + Tpi(s)) + . . . , we find after collecting 
terms of leading order in r that the pair {9o{s), pq(s)) is governed by 



d^n 



ePii9o{s),pois)), 



dpo 



< s < 1, 



as as 
with initial conditions {9q{0), pq{0)) = {9n,Pn)- The solution {9q{1), pq{1)) = {9^,p'^) (obtained 



(21) 



numerically) can then be taken as initial data {9{t),p{t)) = {9^^,p^) in (18)-(19) to form the stroboscopic 
map 



T-jp+{l 



mod 1, 



Pn+1 = P„ e 



(22) 
(23) 



Note that this has been constructed using a matched asymptotic expansion, using (17), and is valid in the 



limit r — !■ 0. For weak forcing, where e <C 1, Pi, 2 vary slowly through the kick and can be approximated by 
their values at (6'„,p„) so that to O(e^) 



0n+l =9n+T + ePi{9n, Pn) + ^ + ^MQn)) (l " e" 

(p„+eP2(e«))e"^^. 



\T\ 



(24) 

Pn+l = {Pn+£i^2[0n}}e - . (25) 

Although this explicit map is convenient for numerical simulations we prefer to work with the full 



stroboscopic map (22)- (23 1, which is particularly useful for comparing and contrasting the behaviour of 



different planar single neuron models with arbitrary kick strength. As an indication of the presence of 
chaos in the dynamics resulting from this system, we evaluate the largest Lyapunov exponent of the map 
( 22][23 ) by numerically evolving a tangent vector and computing its rate of growth (or contraction); see 
e.g. 



32 for details. 



In Fig. 10 we compare the functions Pi 2 for both the FHN and the ML models. We note that P2 for the 



FHN model is near for a large set of 6*, whilst the same is true for Pi for the ML model. This means that 
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kicks in the FHN model will tend to primarily cause phase shifts, whilst the same kicks in the ML model 
will primarily cause shifts in amplitude. 




Figure 10: The blue curves show the change in 6 under the action of a pulsatile kick in u, whilst the red 
dashed curves show the change in p under the same kick. The top plot is for the FHN model, whilst the 
bottom plot is for the ML model. We evaluate the effect of the kicks at p„ = 0, where we observe the largest 
changes in under the action of kicks. 



We plot in the top row of Fig. [TTjthe pair (6I„, 6'„+i), for (|24])-(|25]) for the FHN and ML models. For the 
FHN model, we find that the system has Lyapunov exponent of —0.0515 < 0. For the ML model the 
Lyapunov exponent is 0.6738 > 0. This implies that differences in the functional forms of Pi^2 can help to 
explain the generation of chaos. 

Now that we know the relative contribution of kicks in v to kicks in {0,p), it is also useful to know where 
kicks actually occur in terms of 6 as this will determine the contribution of a train of kicks to the {d,p) 



dynamics. In Figs. 11 C and D we plot the distribution of kicks as a function of 9. For the ML model we 
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A ™N LE=-0.0515 



ML LE = 0.6738 



c 




'i(e,o) 

1.0 



Figure 11: Panel A shows successive iterates of 9 in system (22)-(23) with functions Pi, 2 taken from the 
FHN model, whilst panel B presents the same iterates but for functions Pi 2 from the ML model. Panel C 
shows Pi 2 for the FHN model, where the bold blue line is Pi and the red dashed line is P2- Superimposed 
on this panel is a histogram displaying how kicks are distributed in terms of 9 alone. Panel D shows the 
same information, except this time for forcing functions from the ML model. Parameter values are a — 3.0, 
A = 0.1, £ = 0.1, and T = 2.0. 
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observe that the kicks are distributed over all phases, while for FHN model there is a grouping of kicks 
around the region where P2 is roughly zero. This means that kicks will not be felt as much in the p 
variable, and so trajectories here do not get kicked far from cycle. This helps explain why it is more 
difficult to generate chaotic responses in the FHN model. 

After transients, we observe a 1 : 1 phase-locked state for the FHN model. For a phase-locked state, small 
perturbations will ultimately decay as the perturbed trajectories also end up at the phase-locked state after 
some transient behaviour. This results in a negative largest Lyapunov exponent of —0.0515. We note the 
sharply peaked distribution of kick phases, which is to be expected for discrete-time systems possessing a 
negative largest Lyapunov exponent, since such systems tend to have sinks in this case. The phase- locked 
state here occurs where P2 is small, suggesting that trajectories stay close to the limit cycle. Since kicks do 
not move move trajectories away from cycle, there is no possibility of folding, and hence no chaotic 
behaviour. For the ML model, we observe chaotic dynamics around a strange attractor, where small 
perturbations can grow, leading to a positive largest Lyapunov exponent of 0.6738. This time, the kicks are 
distributed fairly uniformly across 9, and so, some kicks will take trajectories away from the limit cycle, 
thus leading to shear-induced folding and chaotic behaviour. 

4 Discussion 

In this paper, we have used the notion of a moving orthonormal coordinate system around a limit cycle to 
study dynamics in a neighbourhood around it. This phase-amplitude coordinate system can be constructed 
for any given ODE system supporting a limit cycle. A clear advantage of the transformed description over 
the original one is that it allows us to gain insight into the effect of time dependent perturbations, using 
the notion of shear, as we have illustrated by performing case studies of popular neural models, in two and 
higher dimensions. Whilst this coordinate transformation does not result in any reduction in 
dimensionality in the system, as is the case with classical phase reduction techniques, it opens up avenues 
for moving away from the weak coupling limit, where e — ?■ 0. Importantly, it emphasises the role of the two 
functions Pi{0,p) and P2{0) that provide more information about inputs to the system than the iPRC 
alone. It has been demonstrated that moderately small perturbations can exert remarkable influence on 
dynamics in the presence of other invariant structures [IG , which can not be captured by a phase only 
description. In addition, small perturbations can accumulate if the timescale of the perturbation is shorter 
than the timescale of attraction back to the limit cycle. This should be given particular consideration in 
the analysis of neural systems, where oscillators may be connected to thousands of other units, so that 
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small inputs can quickly accumulate. 

One natural extension of this work is to move beyond the theory of weakly coupled oscillators to develop a 
framework for describing neural systems as networks of phase-amplitude units. This has previously been 
considered for the case of weakly coupled weakly dissipative networks of nonlinear planar oscillators 



(modelled by small dissipative perturbations of a Hamiltonian oscillator) [33 - 35 .It would be interesting 
to develop these ideas and obtain network descriptions of the following type 



= 1 + /i {Si , Pi) + w»j iJi {01 , Oj ,Pi,Pj) 
j 

Pi = A{9i)pi + WijH2 {Oj , dj ,Pi,Pj), 



(26) 
(27) 



with an appropriate identification of the interaction functions Hi_2 in terms of the biological interaction 
between neurons and the singe neuron functions Pi. 2. Such phase-amplitude network models are ideally 
suited to describing the behaviour of the mean-field signal in networks of strongly gap junction coupled ML 
neurons 



36 37 , which is known to vary because individual neurons make transitions between cycles of 



different amplitudes. Moreover, in the same network weakly coupled oscillator theory fails to explain how 
the synchronous state can stabilise with increasing coupling strength (predicting that it is always 
unstable), as observed numerically. All of the above are topics of ongoing research and will be reported 
upon elsewhere. 



Appendices 

A Derivation of the transformed dynamical system 

Starting from 

i = f{x) -f eg{x,t), 

we make the transformation x{t) = u{9{t)) + ({9{t))p{t), giving 

'du{d) dc{e) 



-P 



+ ao)p = fH0) + cie)p) + eg{u{e + c(0)p, t). 



de de 

We proceed by projecting (29) onto £,{9), using The left hand side of ([29| now reads: 



du 
d6 



^ deP 



de 
dt' 



where ^'^ denotes the transpose of ^ and the right hand side of ( 29 ) becomes 



ef{u + Cp) + .9(" + Cp) 
du 



de 



^ de" 



+ ef{u + Cp) - e.f{u) - e^^p + ^'^giu + Cp, t). 



(28) 



(29) 



(30) 



(31) 
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Thus, 



1 + /i (0, p) + eh^{0, p)g{u{e) + C(0, t)p, t), 



(32) 



where 



and 



h{e,p) 



Au 



h{e,p) = -h^{e,p)^p(0) + h^{e,p) [fiu + Cp) - /H] . 



Upon projecting both sides of ( [29| ) onto ({0), the left hand side reads 



du dC 



dp 
di' 



(33) 



(34) 



(35) 



whilst the right hand side becomes 

C^fiu + Cp) + <^9{u + CP, t) = -Cfiu) + C^D/C - C^D/C + ef{u + (p) + eg{u + (p, t), (36) 

since C,^ f{u) ~ C,'^du/d6 = and where D/ denotes the Jacobian of /. Putting together the previous two 
equations yields 

dm 



where 



p = Aie)p + f2{9,p) + eC I 



A{e) = e 



de 



de 



-ph{e,p) 



D/c 



g{u{e) + ao)p.t), 



(37) 



(38) 



h{e,p) = -C^^ph + C [f{u + Cp) - f{u) - d/Cp] . 



(39) 



It may be easily seen that fi{9,p) — 0{p) as p — > and that f2{&,0) ~ and df2{d,0)/dp — 0. Overall, 



combining (32 1 and (37) we arrive at the transformed system: 



= 1 + ,h{e, p) + eh^{e, p)g{u{9) + ae, t)p, t), 



p = A{e)p + h{e,p) + ee 



g{u{e) + mp.t)- 



(40) 



In order to evaluate the functions /i, /2 and A for models with dimension larger than two, we need to 
calculate dC,/d9. Defining by ^i{0), the direction angles of ^(6*), we have that 



C. = e. - , (ei + m) = e. - , ' ^^^,L, (d + m) , ^ = 2, 3 



1 + cos 7i 



(41) 
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where the index i denotes the column entry of ( and x ■ y denotes the dot product between vectors x and y. 
Defining 

and 

«^.W = l + eij-?,W, (43) 

where j denotes the row index, we have 

de de ""'de- ^ ' 

By the quotient rule for vectors we find that, 

du, (l + ei)(ei-i)-(e,-$W)(ei§) 



and that 

de de ' 



(45) 



(46) 



Overall, we have that 



da,_ e,-m , / (l + ci)(e..i)-(e.-^(g))(eigg^ 

-df - -i+e,.m ~ [ oT^Tcw? ' • ^ 

B Gallery of models 
B.l Morris-Lecar 

The ML equations describe the interaction of membrane voltage with just two ionic currents: Ca^+ and 
. Membrane ion channels are selective for specific types of ions; their dynamics arc modelled here by 
the gating variable w and the auxiliary functions Woo, Tw, and TOoo- The latter have the form 

moo{v) = ^[l+tanh((t;- wi)/t;2)], Tyj{v) = 1/ cosh ((t; - U3)/(2t;4)) , 

w^{v) = ^[l + ta.nh{{v-V3)/v4)]. (48) 

The function moo{v) models the action of fast voltage-gated calcium ion channels; Vca. is the reversal (bias) 
potential for the calcium current and gca the corresponding conductance. The functions Tyj{v) and w^{v) 
similarly describe the dynamics of slower-acting potassium channels, with its own reversal potential v-^, and 
conductance ^k- The constants wieak and gieak characterize the leakage current that is present even when 
the neuron is in a quiescent state. Parameter values are C = 20.0/(iF/cm^, gi = 2.0mmho/cm^, g^^ = 
8.0mmho/cm2, g^^ = 4.0mmho/cm2, (j) = 0.23, I = 39.5/iA/cm2, vi = -60.0 mV, V], = -84.0 mV, v^^ = 
120.0 mV, vi = -1.2 mV, V2 = 18.0 mV, ^3 = 12.0 mV, and V4, = 17.4 mV. 
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B.2 Reduced Connor-Stevens Model 

For the reduced CS model, we start with the full Hodgkin-Huxley model, with to, n, h as gating variables 
and use the method of equivalent potentials as treated in [26, , giving rise to the following form for the 
function g 



G{v,u) 



dF 
'dh 



hoo{v) ~ hoo{u) 

Th{v) 



dF 
dn 



/ df dhoo{u) ^ df dnoo(u) 
\dhoo du duoo du 



(49) 



where dF/dh and dF/dn are evaluated at h — hrx,{u) and n = n^{u). For the gating variables (a, b), we 
have 



aoo{v) = 
boo{v) = 



1 + e 28.93 / 



1 

" 1. + 53.3 

1 -j- C 14.54 



Ta(v) = 0.3632 + 



1.158 



11 + 55.96 I 
1 + e 20.12 



n{v) = 1.24 + 



2.678 

~ f + 5Q 

1 + e 16.027 



(50) 
(51) 



Parameter values are C = 1 /iF/cm'^, ~ 0.3mmho/cm^, = 36.0 mmho/cm^, = 

47.7mmho/cm2, / = 35.0^A/cm2, = 80.0mV, = -75.0mV, = -77.0mV, = -54.4mV, and 

UNa = 50.0 mV. 



B.3 FitzHugh-Nagumo Model 

The FHN model is a phenomenological model of spike generation, comprising of 2 variables. The first, 
represents the membrane potential and includes a cubic nonlinearity, whilst the second variable is a gating 
variable, similar to w in the ML model, which may be thought of as a recovery variable. The system is 

[IV = v(a — v){v — 1)+/ — w, w = V — bw, (52) 

where we use the following parameter values: fi = 0.05, a = 0.9, / = 1.1, and b — 0.5. 
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List of abbreviations 

• ML - Morris-Lccar 

• FHN - FitzHugh-Nagumo 

• CS - Connor-Stevens 

• LE - Lyapunov exponent 

Additional files 

The additional files for this paper may be found on the additional files page at 

http: / / www.maths.nottingham.ac.uk/personal/pmxkw2 /MultiplePages/Welcome.html 
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